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We calculate the magnon modes in the presence of a vor- 
tex on a circular system, combining analytical calculations 
in the continuum limit with a numerical diagonalization of 
the discrete system. The magnon modes are expressed by 
the ^-matrix for magnon-vortex scattering, as a function of 
the parameters and the size of the system and for different 
boundary conditions. Certain quasi-local translational modes 
are identified with the frequencies which appear in the trajec- 
tory X{t) of the vortex center in recent Molecular Dynamics 
simulations of the full many-spin model. Using these quasi- 
local modes we calculate the two parameters of a 3'^'*-order 
equation of motion for X{t). This equation was recently de- 
rived by a collective variable theory and describes very well 
the trajectories observed in the simulations. Both parame- 
ters, the vortex mass and the factor in front of X, depend 
strongly on the boundary conditions. 

75.10.Hk,75.40.Gb,75.40.Mg,02.60.Cb 



I. INTRODUCTION 

Two-dimensional (2D) magnets have been intensively 
investigated due to different reasons: There are several 
classes of quasi-2D materials for which the magnetic in- 
teractions within planes of magnetic ions are typically 
three to six orders of magnitude larger than the interac- 
tions between the planes. These classes include layered 
magnets (e.g. Rb2CrCl4), graphite intercalated com- 
pounds (e.g. C0GI2), magnetic lipid layers (e.g. man- 
ganese stearate), and high- Tc superconductors. For theo- 
reticians 2D magnets with XY- or easy-plane symmetry 
are particularly interesting due to the existence of vor- 
tices, which have nontrivial properties from the viewpoint 
of certain homotopical groups, see reviews 

There are two types of vortices: In-plane (IP) ones 
with all spins lying in the easy plane, and out-of-plane 
(OP) ones having non-zero spin components orthogonal 
to the easy plane. Both types have a TTi-topological 
charge or vorticity q — ±1, ±2, . . ., which determines the 
directions of the spins in the easy plane far away from 
the vortex center. The OP-vortices have an additional 
7r2-topological charge (Pontryagin invariant, see Ref. ||) 
Q = —\qp where p is integer, p is denoted polariza- 
tion because its sign determines to which side of the easy 
plane the spins point in the vortex center. Therefore 



the IP-vortices can be considered as having p — Q. For 
Q 7^ there is a gyrocoupling force, or gyeoforce, which 
is formally equivalent to the Lorentz force.Q 

Besides the vortices, which are strongly nonlinear ex- 
citations, there are also magnons. The dynamic prop- 
erties of 2D easy-plane magnets can be described by a 
phenomenological theory whick assumes two ideal gases 
for the vortices and magnons.u The vortex-magnon in- 
teraction naturally is the next question which arises. In 
particular, one would like to know how the dynamics of 
the vortices is affected by this interaction. 

For ID magnets the soliton-magnon interaction is 
nearly reflectionless. In this case the main effect of the 
soliton-magnon interaction is the change of the magnon 
density of states, which affects strongly the soliton den- 
sity in thermal equilibrium.cl For 2D magnets the vortex 
density can be obtained from thftjoorrelation length by a 
renormalization group approach,LrLl therefore the vortex- 
magnon interaction is not so important for the density. 
In principle, it can be important for the damping force 
which acts oh a vortex (or soliton) in a near-equilibrium 
magnon gas.B On the other hand, the damping force for 
the 2D case can be obtained by general hydrodjjnamic 
theories for the relaxation processes in magnets.cJ 

The most important effect of the vortex-magnon inter- 
action seems to be that for a finite system certain magnon 
modes are excited due to the vortex motion, and_i£ipi 
versa. Such modes were obtained in recent paperstirE^I 
by a numerical diagonalization for relatively small, cir- 
cular, discrete systems. A calculation in the contifc 
uum limit was presented for the antiferromagnetic case.113 
Analytical investigations were done foij-4)lanar vortices 
in antiferromagnetsta and ferromagnets. Il3 These articles 
demonstrate nontrivial proper|ties of the eigenrnodes, e.g., 
the presence of quasi-locaOE3 or truly localEj modes. 
Moreover, the relevance of these modes for the vortex 
dynamics wa s ah o wn, in particular an IP-OP transition 
was predictedll2lll3 and an effective vortex niass was de- 
fined and calculated using the above modes.Ej However, 
for OP-vortices all this was based on the numerical diago- 
nalization of small systems with fixed (Dirichlct) bound- 
ary conditions (b.c). 

This article deals with a more general theory of vortex- 
magnon coupling for arbitrarily large systems with cir- 
cular symmetry and general boundary conditions. The 
method consists of a natural combination of analytical 
calculations in the continuum limit with numerical diag- 
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onalization of discrete systems. The main issue is the 
calculatioH of the scattering matrix for vortex-magnon 
colhsions.Ol Using these data general formulas for the 
eigenfrequencies are obtained, as a function of the param- 
eters and size of the system, and for different boundary 
conditions. 

An important, point is a link with a recent collective 
variable theoryEZI in which an equation of motion was de- 
rived which describes very well the vortex trajectories ob- 
served in computer simulations. This equation is 3''^-Q&j 
der in time, and is a generalization of the Thiele Eq£3't2l 
which is first order in time. Using the scattering data 
the parameters of the 3''''-order equation of motion can 
be calculated, and, agree very well with those obtained 
from simulations O 



II. THE MODEL AND ELEMENTARY 
EXCITATIONS 

We consider the classical 2D model of a Heisenberg 
ferromagnet (FM) with the Hamiltonian 



H = 



(1) 



{n,n') 



Here J > is the exchange integral, and < A < 1 de- 
scribes easy-plane anisotropy. The spins S are classical 
vectors on a square lattice with the lattice constant aq. 
n, n' denote nearest-neighbor lattice sites. Our main in- 
terest lies in the small anisotropy case which corresponds 
to 1 - A < 1. 

A continuum model for FMs can be derived from (||) 
in the usual way defining the unit vector of magnetiza- 
tion as a function of continuous variables r and i.e. 
rh{f,t) = Sfi{t)/S. The dynamical equations for m have 
the form of the well-known Landau-Lifshitz equation, 
see Refs. In usual angular variables [rux + iruy = 
sin exp (?(/)), TTiz = cos 9] they can be written as 



V'0 



sm tl cos ( 



' V 

V(sin2 eV(j)) 



smt 



cry dt 
sm9de 



cry 



dt 



(2) 
(3) 



where r„ and c are the characteristic length scale and 
magnon speed, respectively. For the Hamiltonian (||), we 
have 



2 



A 



1-A 



2JSaoVl ~ A 



(4) 



where we set h — 1. Note that the Eqs. arise in 

the long- wave approximation (ao|Vm| <C 1) not only for 
the model we are considering here, but for a whole set 
of discrete models, for example, FMs on different kinds 
of lattices and FMs with additional single-ion anisotropy. 
Merely the expressions for c and change. 



For the homogeneous ground state (all spins are par- 
allel and confined to the easy plane) the 2D model has 
well-known magnon excitations with the gapless disper- 
sion law 



u^cHi + k^iy/^, 



(5) 



where /c = |fc| and k is the magnon wavevector. 

Other well-known excitations are the OP-vortices, de- 
scribed by the formulas 



9o(r) 



(6) 



where r and x are polar coordinates and fo is an ar- 
bitrary constant. The function 9o is the solution of 
a nonlinear ordinary jdifferential equation which cannot 
be solved analyticallyfl^H Numerical integration, however, 
shows a surprisingly good agreement with the data ob- 
tained from the numerical analysis of the full discrete 
model (0), even for r„ = l.Soo (A = 0.9), Ref. This 
means that at least concerning static properties strong 
inequalities like ao|Vm| <C 1 can be replaced by usual 
ones. As will be seen, the same is correct for some dy- 
namical properties, too. 9o fulfills the boundary condi- 
tions cos 00 — > p for r — > and ""/^ for r — > oo. 
q and p are topological charges of vorticity and polar- 
ization. In the following OP-vortices with q = 1 will be 
considered; for definiteness we set p = 1. 



III. NORMAL MODES ON THE VORTEX: 
CONTINUUM APPROACH 

We consider small deviations from the static vortex 
solutions, i.e. 



(7) 



Substituting this in the equations (||) and and lin- 
earizing in 'd and /i gives the following set of coupled 
partial differential equations: 



2q cos 9o 9/i 
x'^ dx 

2q cos 6*0 d'd 
x^ dx 



'~'dt 
ry d'd 
'~~dt' 



(8) 
(9) 



where x = r/ry and V^, — r^V. The additional fac- 
tor (sin^o)^^ hi was introduced for convenience, be- 
cause it leads to equations that are symmetric in and 
fi with Schrodinger-type differential operators in front. 
The "potentials" ViCx^, V2 (x) have the same form as for 
the antiferromagnetJi3 



x"^ 



008 26*0, V2^ 



[4-1 






cos^ 6*0 - 




x^ 


dx 



(10) 



but the dynamical parts differ strongly. 
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In order to solve (|[ ^) the following ansatz for i? and 
fi is appropriate: 

+00 

'? = E E {[/; + */a]exp(imx + zc^ai) + c.c.} (11) 

n m— — 00 
+00 

A* = X! X! {[5q +*5a]exp(imx + ic^at) + C.C.} (12) 

n m— — 00 

a = (n, m) is a full set of numbers labeling the magnon 
eigenstates. Substituting this ansatz gives two uncoupled 
sets of two equations for the pairs of functions (/', g") on 
the one and (/", g') on the other side. However, the equa- 
tions for (/", g') can be obtained from the corresponding 
equations for {f',g") simply by replacing m —^ —m and 
Lu — > —Lo. Thus, among the four functions /', /", g' and 
g" only two are independent. Therefore, instead of ( pl| , 
|l2| ) we can use the simplified ansatz 

+00 

■& = ^ ^ fa{r) COs{mX + OJat + 5m) (13) 

71 rn— — 00 

-\-OQ 

M = X! 9a{r) fiin{mx + u;at + 5,n) (14) 



n m— — 00 



with arbitrary phases Sm- For the functions / and g (the 
index a is omitted in the following) we then finally obtain 
the following two differential equations: 





f ^ ^ 


9 




x dx 


x"^ 




1 d 


m 


dx^ 


X dx 


x"^ 
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LOT. 
C 

UJ7\ 
C 



V 



2qm cos t 



2qm cos 0Q 



9 (15) 
/■ (16) 



f and g cannot be determined analytically from Eqs. ( [T^ , 
[iq ). The "potentials" 14(2;) and V2(a;) in these equa- 
tions are not small, and the use of the Born approxi- 
mation looks inadequate; for IP-vortices and m = 0, ±1 
this was already mentioned in Ref. |l^. Only the asymp- 
totic behavior can be calculated. For r — > we obtain 
5j / ~ rl*+'"l , which describes the presence of a "hole" in 
the functions z?, /i at the vortex core for large values of m 
(see Section IV). Only the case |to| = 1 can be considered 
in a long-wave approximation, see below (122 



For the IP-vortex, where cos^o is zero, the modes cor- 
responding to ±|m| are degenerate. As a change of sign 
in the number m can also be interpretated as a change of 
the sense of rotation of the cigenmode (change of sign in 
the eigenfrequency lj), physically we have the situation of 
two independent oscillators rotating clockwise and coun- 
terclockwise with the same frequency (which can also be 
combined to give two linear oscillators in independent 
directions). 

For the OP-vortex, however, with cos^o 7^ 0, this de- 
generacy is removed due to the presence of the term with 
2qm cos 9q jr^ . Therefore the two senses of rotations cor- 
responding to positive and negative frequencies are not 



equivalent anymore. The presence of this term is also 
responsible for the fact that the Eqs. (|l^) and (|l6| ) are 
only invariant under the combined conjugation of both 
topological charges (<7,p) — > {—q,—p). The product qp is 
proportional to the magnitude of the so-called gyrovec- 
tor, which acts on the vortex like a (self-induced) mag- 
netic field in z-direction for a charged particle. There- 
fore, physically, for \m\ > 1 the removal of the ±|m|- 
degeneracy in the OP-case can be understood as the ef- 
fect of a gyroscopic force. 

In order to describe the magnon scattering due to 
the vortex, we note that without the vortex the eigen- 
value problem (EVP) (|, |) can be reduced to a usual 
Schrodinger EVP with the general solution 



^J■ = ^ CmJmikr) sm{mx + tot + dm) 



(17) 



'd ^^Cm—^^===Jmikr)cos{mx + u}t + Sm)- (18) 
V 1 + k^rl 

The Jm are Bessel functions, k and uj are connected by 
the dispersion law (|^) , and Cm and 6m are arbitrary con- 
stants. For a finite circular system with radius L and 
general boundary conditions 

f d'd\ 
a^ + 6r„— j =0 , ("^ + ^''"'^) 

(19) 

with arbitrary constants a and b, the values of the fre- 
quencies d! are determined by the roots of linear com- 
binations of Bessel functions and its derivatives. Note 
that all frequencies for the system without a vortex have 
a 1/L-dependence with respect to the system size. For 
free boundary conditions, i.e. a = in (p^, there is a 
Goldstone mode (GM) with oj ~ due to the rotational 
symmetry in spin space of the model. For a ^ this 
symmetry is obviously broken. 

In the presence of a vortex, one can use the approxi- 
mate formula / « kryg/{l + k^r^)^^'^ for large distances 
r ^ ry and then the EVP (|l^, |l6|) reduces to the form of 
a SEVP, too. The solutions have the same form as ( p^ , 
[1^ ), one only has to replace 



J„i{kr) -> J„i(kr) + a„i{k)Ym^{kr). 



(20) 



Here the Ym, are Neumann functions. The quantity (Tm{k) 
determines the intensity of magnon scattering due to the 
presence of the vortex. In usual scattering theory for 
2D Schrodinger equations the S'-Matrix Sm{k) can be 
represented as 



Sm{k) — 



1 - i<Jm{k) 
1 + i<Jm{k) ' 



(21) 



The values of Umik) are determined by the shape of 
the solution near the vortex core, practically at r/vy < 
3 ... 4. Due to the structure of the continuum equa- 
tions (15, |l^), <7m{k) must be some universal function 
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oi K = kr^, independent of the concrete values of the 
magnetic coupUng constants, the system size or the ap- 
plied boundary conditions. 

o'm(fc) could be calculated directly by solving the set 
of equations ([l5|, |l6|) using a shooting procedure as de- 
scribed in Ref. |13| for the case of the antiferromagnet. 
The shooting parameter is the free constant which ap- 
pears in the general solution of the set of two coupled 
Schrodinger-like EVP. In this article we used another ap- 
proach and extracted the scattering data from the eigen- 
frcqucncies uji we found numerically for the discrete 2D 
system in Section [V. 

In the next section the general properties of am (k) will 
be established for different m, and characteristic features 
of the normal modes will be discussed. Here we consider 
only the case |m| — 1, for which the S-matrix in the 
long-wave approximation can be calculated analytically. 

This can be done by using the fact that for fc = = 
0) and L — > oo a non-trivial solution of (^5|, |l^) is known, 
and by applying a special perturbation theory. This solu- 
tion reads /o = dOg/dx, go = — ainOo/x, and corresponds 
to a vortex displacement (translational Goldstone mode) . 

In order to construct the asymptotics of such a solution 
for a small but finite frequency, we make the ansatz 

/ = ideo/dx)[l + a{x)] , g = ~{sm9„/x)[l + (22) 

where a{x),P{x) must be proportional to lo. 

Inserting this ansatz into the set of Eqs. (jlj, mul- 
tiplying jig ) with go and ([l5| ) with /q, and adding the 
results, one obtains the equation 



sin ^0 



2^(cos0o)'. (23) 
c 



Here the prime denotes d/dx, and small terms like 
wa, uj(3 were omitted. The formal solution of this equa- 
tion can be written as 

(3{x) = /3(0) - / a'{x){xe'Q/ sinOofdx 
Jo 

-{ojrjc) [ [x/cos^{0Q/2)]dx. (24) 
Jq 

It is easy to see that the last term is divergent for x ^ 0, 
while the integral with a'{x) is convergent due to the 
presence of 0'{x) ~ exp (— x). Thus, far from the vortex 
core we simply have 

l3{x) ~ — ((jjrt,/c)a;^ = —kr^x'^, (25) 

valid in the region 1 <C x <C (l/fcr^,)^/^. g{r) then reads 



g{r) ~ — ~ kr^ ( — 
r \ r„ 



(26) 



On the other hand, in the region r (ry/ky^'^ 

the arguments of the Bessel functions in (^) are small: 



kr ^ (fcri,)-^/^ <^ 1. Using the asymptotics Ji(fcr) = 
kr /2,Yi(kr) — —2/nkr, one can rewrite Eq. (^ ) in the 
form 



5(0 



kr nkr. 

const < — -\ — , , 

2 4 V Trfcr 



Ji(fcr) 



irkr,, 



'-Yi{kr), 



(27) 



which gives us the long-wave approximation of the scat- 
tering amplitude: cTipi dznkry/A. This is in good 
agreement with the numerical data as we will see in the 
next section. 



IV. NUMERICAL DIAGONALIZATION AND 
SCATTERING DATA CALCULATION 

A numerical diagonalization of the same model we are 
considerina-iere was also already performed by Wysin 
and Volkel.td These authors found, for example, the split- 
ting of the TO 7^ modes for OP vortices, as discussed 
above, and the lack of such a splitting for IP vortices. 
For a detailed description of the eigenvalues and eigen- 
functions of (|l|) we therefore refer to their paper. For the 
rest of this section we relate some points which were not 
discussed in Ref. |l^ including a short description of dis- 
creteness effects. Our main interest however, lies in the 
scattering of magnons due to the presence of a vortex. 

We consider a circular system of radius L on a square 
lattice with either fixed (Dirichlet) or free (Neumann) 
boundary conditions. Static vortex solutions were ob- 
tained using a relaxation procedure similar to that de- 
scribed in Ref. ^ or, in slightly more detail, in Ref. 

One problem with the diagonalization of 2D systems 
is that the number of sites in the system and therefore 
the dimension of the matrix to diagonalize grows very 
rapidly with the linear size of the system. Already for a 
rather small system with radius L — 20ao this matrix has 
dimension 2528 x 2528. If we use a conventional numeri- 
cal method which needs to store the whole matrix in the 
memory of the computer, this is already more or less the 
size limit for matrices which can practically be diagonal- 
ized. In order to go further, we applied a special method 
which takes advantage of the structure and sparseness of 
the matrix. Details are described in Appendix A. Us- 
ing this method we calculated the lowest eigenvalues and 
eigenvectors for systems up to L/a^ = 100 corresponding 
to a matrix with dimensions 62856 x 62856. 

The eigenvalue spectrum we obtained for fixed bound- 
ary conditions is plotted in Fig. |], as a function of the 
anisotropy parameter A, for different values of the an- 
gular momentum quantum number |m|. For compari- 
son, we have also calculated the eigenvalues for a system 
linearized around the ferromagnetic ground state of the 
model (i.e. without a vortex being present in the system). 
These data are the dashed lines in Fig. |l|. Ac « 0.70 is 
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the critical value of anisotropy which separates the dif- 
ferent regimes of stability of IP- and OP-vortices. As 
can easily be seen, the mode which is responsible for this 
transition is the m = mode (the lowest branch in Fig. 
|l|), which becomes soft at Ac- Mainly the m = and 
|m| = 1 modes are sensitive to the IP-OP-transition. The 
eigenvalues above Ac are split where for fixed boundary 
conditions the lower branch corresponds to negative m 
and the upper branch to positive m. 

For A < Ac, only eigenmodes with odd |to| are twofold 
degenerate, as predicted in the continuum theory. For 
even |to| the expected degeneracy is slightly removed. 
This is a discreteness effect, due to the fourfold symme- 
try of the underlying square lattice: eigenvectors with 
even |m| tend on the one hand to align with a coordinate 
system parallel to the symmetry axes of the underlying 
lattice, and on the other hand with a coordinate system 
parallel to the lattice diagonals. The latter eigenmodes 
oscillate with a slightly higher frequency. 

The classification of the eigenvalues according to val- 
ues of |m| was done by expanding the eigenvectors in a 
Fourier series similar to (|l^, |lj) (excluding the time de- 
pendent part). For A = 0.9 some of the radial profiles of 
eigenvectors obtained this way are plotted in Fig. ^j. The 
numerical computation always gives a complex eigenvec- 
tor. For odd |m|, the real and imaginary radial profiles 
coincide with a high numerical accuracy. For even \m\ 
however, as a consequence of the discreteness effect al- 
ready described above, the radial eigenfunctions of the 
real and imaginary part differ slightly, especially near 
the center of the vortex. (Fig. ^ only shows a mean value 
in these cases.) Also, for even |m|, the sense of rotation 
of the eigenmode can only approximately be determined: 
modes in the continuum theory which belong to positive 
and negative m intermix in the discrete system for even 
|m|. This second effect is more pronounced for larger \m\ 
and larger system sizes. As a consequence, these eigen- 
modes are not in the separated form of Eqs. (^3|) and 

The main point in our approach is to calculate the 
scattering matrix Sm{k) and to use its values for the in- 
vestigation of general properties of the magnon modes 
in the presence of the vortex. For the extraction of the 
5m (fc) dependence on the basis of Eq. (|l^, ^) we only 
need values of eigenfrequencies for one boundary condi- 
tion for different values of the system size. Concretely, we 
used the fixed b.c. (& = in (|ig|)) leading to the simple 
formula 



(28) 



with ki = kiiuji) according to the dispersion law 

We calculated the scattering data for four different val- 
ues of A, namely A = 0.80, 0.85, 0.90 and 0.98, and for 
different system sizes up to L/oq = 100. These data are 
presented in Fig. ^ as a function of k = fcr„. We note 
that the points am{n) fit well to one curve for different L 



and A, at least for small values of ka^. (Strictly speaking, 
the inequality fcao <C 1 should be satisfied, but even for 
/cao ~ 0.4 the discrepancies among different values of A 
are negligible for the mainly interesting case |m| = 1.) 

For the general boundary conditions ( p^ ) the frequency 
spectrum is determined by the dimensionless equation 



a[J^{kL) ■ 
kr,b[J'^{kL) 



a, 



,(kr,)Y^^{kL)\ 
,{kr,)Y:^{kL)] 



0. 



(29) 



We start our discussion for the modes with |m| > 1 and 
m = 0. The translational modes (|m| = 1) will be con- 
sidered in the next section. 

For |m| > 1 and fixed b.c. (6 = in (|l^)), taking into 
account <t,„ < 0, it can easily be seen that Eq. (|29| ) has 
no solution in the region kL <C 1 where Jm{z) ~ (z/2)™ 
and Ym{z) ^ — (2/z)''"/7r. The n'th root can be written 



as kL — 



with y„ 



< Jr. 



Here y^,™ and 



jm,n are the n'th zeros ofYm{z) and Jm{z), respectively. 
Because jm,n and ym.n are smaller than 7r(n -t- m)/2 for 
large n or m, the low-lying modes have the values fc„^m ~ 
max(n, m)/L. This means am 1 for r^/L <C 1, hence 
the eigenfrequencies can approximately be described by 
the zeros of Jm{z) (which have approximately the same 
values as for the system without vortex). This can also 
be observed very clearly in Fig. |l . 

For the case a ~ 6 in the b.c. (19^ the same holds be- 
cause of the presence of the small parameter kry in the 
term with the derivatives. For a = (Neumann b.c), 
however, the situation changes. First of all, the formula 
for kn^m can be presented in the form fc„^m = j'^.n/^ 
where j'^ „ is the n'th zero of J^iz). More importantly, 
a new root of ( p9| ) appears, with k = k^^n ^/L 
corresponding to the lowest mode with \m\ > 1. In 
order to show this we use again the asymptotics of 
Jm{z) and Ym{z) for z ^ 1. Applying the formula 
Cm = — ^m(fcfu)", one easily obtains 



l^vko rn — 



p/2ni 



2ry 
L 



1 



P = 



1 — a/2m 



> 1. 



(30) 



The corresponding frequency cjq.tjj. ~ cfco,m decreases 
faster than 1/L for L ^ oo. 

The eigenmodes with m — can be described in the 
same way as for |m| > 1. Only the values kn^ — 
{l/L)jo^n or knfi — (l/_L)jQ„ are present in this case, 
and kn,o ~ ^/L for general boundary conditions. The 
only exception is the lowest mode for Neumann b.c. For 
this mode the analysis of the asymptotics of cylindrical 
functions gives exactly a; = 0. This value corresponds to 
the Goldstone mode due to the presence of a rotational 
symmetry in spin space which is still exact for the finite 
system. 
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V. TRANSLATIONAL MODES AND 
COLLECTIVE VARIABLE APPROACH 



renerate an image of opposite vorticity. Therefore from 
j33|) two equations follow: 



This section deals with the analysis of the translational 
modes |m| = 1 and the application of their properties for 
the construction of equations of motion for the vortex 
dynamics. 

From Fig. || one can see that the scattering parameter 
for |m| = 1 have a linear dependence (J±i{k) — ^Aikr^. 
The long- wave approximation (^) gives Ai ~ (7r/4), 
which is in good agreement with Fig. 0. Using the asymp- 
totics Ji ~ z/2, Yi ~ —2/ttz for z <C 1, we get for the 
frequency of the lowest translational mode loq = ztcry/L'^ 
where the positive sign corresponds to Neumann and the 
negative sign to Dirichlet boundary conditions, respec- 
tively. Using cry — a^JSVX and replacing y/X 1, 
which is natural in the continuous small anisotropy ap- 
proximation, we obtain luq = ^JSuq/L^. 

This analytical result is in agreement with the eigen- 
value data obtained numerically. However, in the discrete 
system, uq does actually not tend to zero for L — > oo 
but to a finite constant cop. Using the collective variable 
approach this constant can easily be understood as a pin- 
ning effect as we would like to show now. It is known that 
the dynamics of the OP vortex can be described j n| | l0 |west 
approximation by the so-called Thiele Equatior£3ll3 



G X X F. 



(31) 



Here G is the gyrovector G = {2TTqpS/al) z, where z is 
the unit vector along the hard (z) axis, X the position of 
the vortex center, and F the force acting on the vortex. 
In the case where the vortex is only slightly displaced 
from the center of the system, Thiele's Eq. just describes 
the dynamics of the lowest translational mode. There are 
two types of forces acting on the vortex. The first one 
is due to an image vortex which resides at the position 
{L"^ / X'^)X , X = \X\, and the second one is due to the 
pinning potential of the underlying lattice. Putting both 
together we obtain 



Gx X 



2TrJS^qq.,X 



kpX, 



(32) 



where qi is the vorticity of the image vortex. We have 
assumed that the pinning potential is harmonic in lowest 
approximation leading to a linear force with spring con- 
stant kp. In the limit X ^ 0, Eq. ( p2[ ) has the solution 
Xi{t) — X cosujQt, X2{t) —XsinujQt with 



UJQ = - 



qqiJSal kpU^ 



L2 



27rS 



(33) 



which, apart from the constant part, agrees with the de- 
pendency we have found above. In order to calculate kp 
from the numerical eigenvalue data one can use the fact 
that Dirichlet boundary conditions generate an image of 
the same vorticity while Neumann boundary conditions 



l-^O, Dirichlet 



^0, Neumann ^" 



JSal 
JSal 



kpCiQ 
2^ 
kpOjQ 
2^' 



(34) 



Adding them leads to the following formula for ujp which 
is independent of the concrete size of the system: 



kpCiQ 
2^ 



^Q, Dirichlet ~^ ^0, Neumann 



(35) 



Assuming a pinning potential of the form 
£'p[cos(27ra;/ao) -|- cos(27ry/ao) -l-2]/2, uip is connected to 
the pinning energy Ep by the formula 



Ep — 



S\LUp\ 



(36) 



We have evaluated (|35|) and (36) for some values of A, 
cf. Table VI. Naturally, the pinning energy decreases 



with decreasing anisotropy because the size of the OP- 
structure increases strongly, making discreteness effects 
less and less important. 

The pinning energy can also be measured directly, com- 
paring the energy of static vortices centered on different 
lattice coordinates. These numbers, labeled E: 



direct 



Table VI , compare very well with Ep due to formula (|3£ 

Next, we discuss the higher modes obeying the |m| = 
1 symmetry. Due to the general features discussed in 
the previous section, these modes can be considered as 
doublets, with frequencies a;„ > and tOn+i < with 
|[j„+i| ~ LOn- For general b.c. the mean frequency can be 
written as 



(37) 



where Xn is the n-th root of the equation aJi{x) + 
bx{ry/ L)J[{x) = 0. For fixed and free b.c. one has 
^n — cji.n/L and lU„ — cj'i n/L, respectively. 

We concentrate on the lowest doublet with the frequen- 
cies uji and UJ2 — —loi on the following reasons: (i) these 
modes can be compared with long-time computer simu- 
lations of the vortex motion and (ii) they can be used for 
the calculation of parameters in equations of the vortex 
motion. In the lowest approximation in r„/L one can 
write 



UJ-l , — 

- . 2x,VT 



-r-Oo 



(38) 



J'l.i 



where Xi = ji^i ~ 3.832 for fixed b.c. and xi 
1.8412 for free ones. For general b.c. xi lies between these 
values. Note, however, that for a ^ b the value of xi is 
still close to ji^i, and the "switching" to the value j[ i oc- 
curs only for small a < b{aQ/L). We remind that ZJi does 
not depend on the scattering data in the first approxi- 
mation in ao/L. The use of scattering data becomes im- 
portant however for the calculation of the doublet width 
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Aw = |w2| — t^i- Considering the leading approximation 
kry ~ cio/L, the width can be presented in the form 



Aw 



n^/X uq 2 xKaxi — b)Yi + xibYQ 

2 ^L^ xi{xia - b)Jo + [{2 - xl)b - xia]Ji' 

(39) 



where Jq,Ji,Yq and Yi aU have the argument Xi. 
One can see that Aw is practically independent of the 
anisotropy (we will omit the coefficient a/A below), and 
it is inversely proportional to L^. This is the same 
dependence as for the lowest translational mode ujq — 
JS{ao/L)^, but the coefficients are different. 

Next we compare our results for luq and the lowest dou=. 
blet wi_2 with data from recent computer simulationstll 
for the motion of one OP-vortex on large circular systems 
(up to L = 72ao) with free boundaries. The Landau- 
Lifshitz Eq. was integrated for a square lattice with 
flo = 1 and the Hamiltonian ([l]), where J — S = 1 and 
A = 0.9. The vortex center performs oscillations around 
a mean trajectory X'^{t) which is a circle with radius 
i?o around the circle center (Fig. ^). The mean trajec- 
tory can be interpreted as a stationary solution of the 
Thiele Eq. (^Tj), where the vortex is driven by the inter- 
action with its image vortex. The rotational frequency 
Wo goes to a constant for i?o ^ 0; the extrapolated value 
0.201 10-3 agrees weh with our result 0.193 10-^ (both 
for L = 72). 

The Fourier spectrum of the oscillations around X"(t) 
shows a doublet to'i 2 with about equal amplitudes, and 
phases =F7r/2 between the two components. As the spec- 
tra in Ref . |l^ were evaluated in a moving polar coordinate 
frame we must add and subtract wq in order to compare 
with our results in Eqs. (BSh and (pS) 



Wl,2 



(40) 



The data agree very well, e.g. for L = 72 and Rq — 
16.11 the frequencies wq = 2.05 10^^, uj[ = 1.546 10"^ 
and W2 — 1.661 10"^ were observed which yields w = 
1.604 10"^ for the mean and Aw = 0.74 lO'^ for the 
difference. This agrees with our theoretical values ZJ ~ 
1.617 10-2 and Aw = 0.7657 10"^ within 0.8 and 4%, 
respectively. 

Note that the very good agreement of our normal mode 
approach with the data from the simulations is not a triv- 
ial result: we have calculated the frequencies appearing in 
the vortex dynamics for free boundary conditions, using 
as input the scattering data for fixed boundary conditions. 
This shows that our scattering theory actually works for 
the general case. 

We conclude that the vortex motion is accompanied 
by, or generates, the two quasi-local magnon modes with 
|m| = 1. These show up in two ways: (i) in the trajec- 
tory of the vortex center as oscillations around a mean 
trajectory, (ii) as oscillations of the dynamic parts of the 
vortex structure (see below) . 



There remains the question whether the trajectory 
of the vortex center can be obtained from an equa- 
tion of motion. If the rigid-shape assumption of Thiele 
is dropped by allowing for a deformation due to the 
velocity,! a 2'"*-order equation of motion can be derivedEll 

which exhibits an additional term MX with an effec- 
tive vortex mass M. The additional term causes cy- 
cloidal oscillationsi with frequency Wc = G/M around 
a mean trajectory. For the vortex mass AI ^ InL was 
predicted, which was not confirmed by coixmuter simula- 
tions, though Refs. |2 23, Recently WysinEfl proposed to 
calculate M by using the two lowest eigenfrequencies, wq 
and wi in our notation. His formula A-Iw = G/ (wq -I- wi) 
gave a linear L-dependence when the frequencies from 
his numerical diagonalization were inserted. The same 
dependence is obtained (in 0{aa/ L)) if we insert our an- 
alytical results for the frequencies. 

However, meanwhile new simulations had been per- 
formed which resulted in the observation of the above 
mentioned doublet wi_2, instead of a single frequency Wc. 
This dynamicSpCan be fully described by a 3'''^-order equa- 
tion of motionllZI which was derived by a collective vari- 
able theory, starting from a generalized traveling wave 
ansatz 



S{f,t) = S{r-X,X,X) 



(41) 



where the vortex shape is assumed to depend also on the 
acceleration. In fact, this dependence can be seen in sim- 
ulations with free b.c. when one considers the spin config- 
urations at the turning points of the trajectory where±he 
acceleration is maximum while the velocity is small.Ej 
The 3'"''-order equation has the form 



GzxX + MX -Gx X = F 



(42) 



with G3 = G^Cz and a new parameter G3. We note that 
this is the most general 3''''-order equation for the given 
easy-plane symmetry. In Ref. ^7[ G3 was defined as an 
integral over r which could not be performed because 
the dynamic vortex structure is not known analytically 
for the vortex core. But the size dependence G3 ^ 
was obtained from the outer region of the integral. We 
will see below that our theory allows for the calculation 
of G3 and M. 

We consider Eq. (^) for the case of vortex motion 
close to the center of the circle where the image force is 
approximately linear 



F ^ -27rqq,X/L^ 



(43) 



Then Eq. ( [i^ can be solved by a harmonic ansatz, where 
the frequency w fulfills 



Gsw^ - Mw^ + Gw = -2TTqqi/L^ 



(44) 



The parameters can be obtained from the three roots w^, 
using Vieta's rules 
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UJ0UJ1UJ2 



-M/G3 
-G/G3 
27rgg,/(G3i'). 



(45) 
(46) 
(47) 



We now identify ujq with the Goldstone mode —qqi/L^, 
neglecting the frequency Up from the pinning force since 
Eq. (^) was derived in the continuum hmit; uji and 
u!2 < are identified with the lowest doublet calculated 
in Eqs. (MM. n 
From Eq. ([47]) we obtainH 



with uj = Y^aJija^ — 5(^1 + \oJ2\), which gives 



G, = 



2x2(1- A) 



(48) 



(49) 



The result 4.634^^ for free b.c. agrees very well with 
4.67ilj''°2 for large L obtained from the simulation 
dataO For fixed b.c, G3 ~ 1.07L^ is about four times 
smaller. 

From Eq. (p5|) we obtain 



(50) 



with Ao; ~ \uj2 \ —LJi- Here all L-dependencies just cancel 
in the lowest approximation on L, thus the vortex mass 
in Eq. ( ^ ) is independent of the. system size, as already 
obtained from the simulations .EZl Inserting Eq. ( ^9|) we 
get for free b.c. 



M 



2(1 -A) 



l)Jl 



(51) 



The numerical value M = 14.74 agrees well with M = 15 
which we have extrapolated for i?o ^ from the data in 
Ref. 17 for three different system sizes {L = 24,36,72). 
We note that the vortex mass (51) is in the same order 



of magnitude as the two-dimensional soliton mass M ~ 
EqI(? in Ref. where Eq = 47r (in units JS"^) is the 
Belavin-Polyakov energy. 
For fixed b.c. 

M 



2(1 - A) V2 a;iJo 



vr Fi 



(52) 



yields M = 7.661, which is about one half of the above 
value for free boundary condition. 

Thus in the equation of motion ( ^ ) only the gyrovec- 
tor G is an intrinsic property of the vortex. The quan- 
tities M and G3, which are connected to the quasi- local 
modes with frequencies wi.2, are determined by the whole 
system "vortex plus magnons" which includes the geom- 
etry of the system and the boundary conditions. For G3 
this is obvious because it strongly depends on L; but M 
does not depend on L (in the lowest order), however, it 



depends on the boundary conditions. On the other hand, 
the S-matrix naturally is determined only by the region 
of the vortex core and can thus be used for the calcula- 
tion of the parameters M and G3 for arbitrary geometry 
and boundary conditions. 

Finally we point out that Eq. ( ^2|) belongs to a whoia 
hierarchy of equations of motion which can be derivedllll 
by taking into account higher and higher time derivatives 
of X{t) in the generalized traveling wave ansatz (|4l|). 
Moreover, only the odd-order equations of this hierarchy 
represent valid approximations, because the even-order 

equations have a very weak leading term, e.g. MX in 
the 2"''-order equation. Therefore the solutions of the 
latter equations are qualitatively different from the solu- 
tions of the odd-order equation and are in fact not con- 
firmed by the simulations (cf. the above discussion of the 
single frequency lJc and the doublet 0^1,2, which exhibit 
quite different dependencies on the system size). The 
5*^-order equation of the hierarchy predicts a second dou- 
blet (^3.4 which was in fact observed in the simulations, 
but only for specially designed initial and boundary con- 
ditions-because the amplitudes of this doublet are very 
small.tZl Thus the 3'''^-order equation represents already a 
very good approximation. The observed additional dou- 
blet naturally also appears in the results of our numerical 
diagonalization in Section [V and its frequencies can be 
calculated by ZZ;„ and (Aw)„, see Eq. (^ and below. 



VI. CONCLUSION 

We have developed a general theory which allows cal- 
culation of the magnon modes of a circular easy-plane 
ferromagnet in the presence of an out-of-plane vortex. 
The method consists of a combination of numerical di- 
agonalization of the discrete system with analytical cal- 
culations in the continuum limit. The frequencies of the 
magnon modes can be expressed in terms of the functions 
<7m{kry), which are independent of the magnetic coupling 
constants, the system size and the boundary conditions. 
The am describe the intensity of the magnon scattering 
due to the presence of the vortex. 

The translational modes with |m| = 1 are particularly 
interesting for two reasons: (i) Their frequencies are iden- 
tified in the vortex motion which was observed in simula- 
tions where the Landau-Lifshitz Equation was integrated 
for the circular discrete spin system, (ii) Using these fre- 
quencies, one can calculate the two parameters of a 3'"''- 
order equation for the vortex motion (a generalization of 
the Thiele Eq.), which was derived by a collective vari- 
able theory starting from a generalized traveling wave 
ansatz. Our calculated parameters agree very well with 
those obtained by describing the simulations using the 
3»'<i_order equation of motion. Both parameters, the vor- 
tex mass M and the 3'''^-order gyrocouphng constant G3, 
depend strongly on the boundary conditions. This is due 
to the fact that both the in-plane and the out-of-plane 
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structure of the moving vortex are not localized. In fact, 
the dynamic parts of the vortex structure oscillate with 
the frequencies of the translational modes with \m\ — 1. 
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APPENDIX A: NUMERICAL 
DIAGONALIZATION OF HAMILTONIAN 
MATRICES 

If a classical Hamiltonian system is been linearized 
around a static solution, an eigenvalue problem 



Az = Xz 



(Al) 



results with a coefficient matrix A fulfilling the general- 
ized symmetry relation 



J is the "symplectic unit matrix" , i.e. 



1 
1 



(A2) 



(A3) 



For obvious reasons the matrix A is called a "Hamilto- 
nian matrix" . It is easy to see that A can always be 
written as JH where H is a symmetric matrix. Note 
that the eigenvalues and eigenvectors of A are complex 
in general. 

For the most general case, Hamiltonian matrices have 
no special properties which could profitably be used in 
the numerical diagonalization. An exception, however, is 
the case of a positive definite matrix H, for which ( |Al[ ) is 
equivalent to an Hermitian eigenvalue problem with pure 
imaginary eigenvalues. In order to show this, we note 
that a positive definite matrix can always be written as 



H = LL' 



(A4) 



the so-called Cholesky decomposition of H, cf. Ref. 
Substituting this in (Al) and additionally defining A :~ 
iui we obtain 



(zLML)(L^z) = -t^(L^z) 



(A5) 



Now, J is antisymmetric and therefore iL'^JL is Hermi- 
tian. 

Unfortunately, it is somehow difficult to use this equiv- 
alence directly for the purpose of numerical diagonaliza- 
tion. The Cholesky decomposition of a sparse symmetric 



matrix can be computed easily only for certain cases, 
otherwise it is very time consuming. For this reason we 
did not directly use formula (A4). Instead we utilized 
a method called Wielandt's version of the inverse itera- 
tion procedure. The basic strategy is to multiply a (ran- 
domly chosen) vector over and over by the inverse of the 
spectral-shifted matrix A, i.e. (A — AI)~^. The resulting 
series of vectors converges to an eigenvector of A, usu- 
ally the one corresponding to the eigenvalue closest to the 
chosen spectral shift A. The details of the method can be 
found in Ref.Ej. For Hamiltonian matrices, supposing a 
positive definite matrix H and therefore imaginary eigen- 
values and an imaginary spectral shift A = jw, Wielandt's 
inverse iteration amounts to the following iteration for- 
mula: 



H 

u}3 



-ujJ 
H 



(A6) 



x^^^ and y^^^ are the real and imaginary part of the 
eigenvector in the j'th iteration step. The initial vec- 
tor {x'--^^ , yf"' ) must be choosen randomly. The matrix V 
is a symplectic projection operator which is defined as 



k<j 



f(fc)Tjy'(fc) 



-J. 



(A7) 



The sum runs over all eigenvectors which were computed 
in previous runs. The purpose of this operator is to avoid 
that the method converges to an already-known eigenvec- 
tor. 

After some iteration steps the parameter uj can also be 
iterated, for example according to the formula 



1 z+Az 



z = x^-'^ 



(A8) 



As z converges to an eigenvector of A, A converges to the 
corresponding eigenvalue. 

A positive definite matrix H is not a principal pre- 
sumption for the inverse iteration, i.e. the above itera- 
tion formula can in principle be generalized to include 
the case of non-definite matrices H. How ever, first of 
all, with H, the coefficient matrix in (A6) is also pos- 
itive definite. Th is a llows to use an efficient numerical 
method to solve (A6). To be precise we have used the 



NAG library function F04MBF for this purpose which is 
based on a Lanzcos algorithm. The NAG routine makes 
it only necessary to supply a "matrix x vector" function. 
Therefore one has to store only the elements of A which 
are non-zero which precisely makes the method suitable 
for large and sparse matrix equation. A second point is 
that for non-definite H rather serious problems with the 
numerical stability of the inverse iteration arise. 



9 



* Permanent Address: Institute for Metal Physics, Ukraine 
Academy of Sciences, 36, Vernadskii St., Kiev 252142, 
Ukraine. 

^ For antiferromssnets this matrix was calculated by a dif- 
ferent method. c3 

* A mass, due to a velocity induced shape deformation, 
was already pjtroduced for the solitons in 2D easy-axis 
ferromagnets 

^ These trajectories correspond to the cycloids with fre- 
quency eB/{Mc) of an electric charge e in a plane with 
a perpendicular magnetic field B and an in-plane electric 
field E, where c is the speed of light. 

^ Eq. ( |46[ ) gives the same result, because a term loqAlo/lu^ 
can be neglected. 
^ A. M. Kosevich, B. A. Ivanov, and A. S. Kovalev, Phys. 
Rep. 194, 117 (1990) 

V. G. Bar'yakhtar and B. A. Ivanov, in Soviet Scientific 
Reviews, Section A. Physics, I.M. Khalatnikov (ed.). Vol. 
16, No. 3 (1993) 
^ B. A. Ivanov and A. K. Kolezhuk, Low Temp. Phys. 21, 
275 (1995) 

* A. P. Malozemoff and J. C. Slonczewski, Magnetic Domain 
Walls in Bubble Materials, Applied Solid State Science, 
Supplement I, Academic Press, New York (1979) 

^ F. G. Mertens, A. R. Bishop, G. M. Wysin, and C. Kawa- 
bata, Phys. Rev. Lett. 59, 117 (1987); Phys. Rev. B39, 591 
(1989) 

•^H.-J. Mikeska and M. Steiner, Adv. Phys. 40, 191 (1991) 
J. M. Kosterlitz and D. J. Thouless, J. Phys. C6, 1181 
(1973) 

* J. M. Kosterlitz, J. Phys. CT, 1046, (1974) 

" V. L. Berezinskii, Sov. Phys. JETP 32, 493 (1970); 34, 610 
(1972) 

E. G. Galkina and B. A. Ivanov, JETP 82, 1168 (1996) 
" G. M. Wysin, Phys. Rev. B49, 8780 (1994) 
^2 G.M. Wysin and A. R. Volkel, Phys. Rev. B52 , 7412 

(1995); B54, 12921 (1996) 
" B. A. Ivanov, A. K. Kolezhuk, and G. M. Wysin, Phys. 

Rev. Lett. 76, 511 (1996) 

B. V. Costa, M. E. Gouvea, A. S. T. Pires, Physics Letters 
A156, 179 (1992) 

A. R. Pereira, M. E. Gouvea, A. S. T. Pires, Phys. Rev. 
B54, 6084 (1996) 

G. M. Wysin, Phys. Rev. B54, 15156 (1996) 

^'^ F. G. Mertens, H. J. Schnitzer, and A. R. Bishop, Phys. 
Rev. B56, No. 5 (1997) 

A. A. Thiele, Phys. Rev. 30, 230 (1973); J. App. Phys. 45, 
377 (1974) 

D. L. Huber, Phys. Rev. B26, 3758 (1982) 

H. J. Schnitzer, Ph. D. thesis (in German), University of 
Bayreuth (1996) 

G. M. Wysin and F. G. Mertens, in Nonlinear Coherent 
Structures in Physics and Biology, M. Remoissenet and M. 
Peyrard (eds.). Springer Lecture Notes, Vol. 393, Springer, 
Berlin (1991) 

G. M. Wysin, F. G. Mertens, A. R. Volkel, and A. R. 
Bishop, in Nonlinear Coherent Structures in Physics and 
Biology, K. H. Spatschek and F. G. Mertens (eds.). Plenum 
Press (1994) 

22 F. G. Mertens, G. M. Wysin, A. R. Volkel, A. R. Bishop, 



and H. J. Schnitzer, in Nonlinear Coherent Structures in 

Physics and Biology, edited by K. H. Spatschek and F. G. 

Mertens, Plenum Press (1994) 
23 A. R. Volkel, G. M. Wysin, F. G. Mertens, A. R. Bishop, 

and H. J. Schnitzer, Phys. Rev. B50, 12711 (1994) 
2** B. A. Ivanov and V. A. Stephanovich, Phys. Lett. 141A, 

89 (1989) 

2^ J. Stoer, and R. Bulirsch, Introduction to numerical anal- 
ysis, Springer, Berlin (1993). 



TABLE I. Dependence of pinning frequency and energy on 
the anisotropy parameter A. 



A 


M/JS 


Ep/{JS^) 


^^'■""V(js2) 


0.80 


0.0672 


0.0214 


0.0233 


0.85 


0.0181 


5.76 10"^ 


5.87 10"^ 


0.90 


1.88 10"=^ 


5.98 10"'' 


6.05 10"* 


0.98 


8.07 10"^ 


2.57 10"^ 


3.85 10-** 
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FIG. 1. Eigenvalue spectrum as a function of the 
anisotropy parameter A and the "quantum number of angu- 
lar momentum" m. Calculated for a circular system of radius 
L = 20 with fixed boundary conditions. The dashed lines are 
eigenvalues with no vortex in the system. 
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FIG. 2. Radial part of eigenfunctions for the lowest eigen- 
values at A = 0.9 using fixed boundary conditions. Solid and 
dashed lines correspond to eigenfunctions with positive m and 
negative m, respectively. 
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FIG. 3. Scattering data for different values of A (□ = 0.80, 
A = 0.85,+ = 0.90, O = 0.98) and different system sizes in 
the range L = 15 . . . 100. 
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FIG. 4. Trajectory of a vortex obtained by numerically in- 
tegrating the Landau Lifshitz equation in time. The simula- 
tion was performed on a circular system with radius L = 36oo 
and free boundary conditions. 
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